Fix classification regression and reproduce JPP 2020 paper results#323
Open
krystophny wants to merge 13 commits intomainfrom
Open
Fix classification regression and reproduce JPP 2020 paper results#323krystophny wants to merge 13 commits intomainfrom
krystophny wants to merge 13 commits intomainfrom
Conversation
ⓘ You are approaching your monthly quota for Qodo. Upgrade your plan Review Summary by QodoRestore volume sampling and bminmax output for classification
WalkthroughsDescription• Restore volume sampling for num_surf=0 to enable classification diagnostics • Fix fast_class and tcut compatibility by removing incorrect mutual exclusivity check • Generate bminmax.dat output file with 101 rows for classification plotting • Migrate classification binning from Fortran postprocessor to Python script • Fix hardcoded array index and VMEC file reference in classification example Diagramflowchart LR
A["startmode=1<br/>num_surf=0"] -->|"volume sampling<br/>uniform s in [0,1]"| B["sample_particles"]
B --> C["compute_pitch_angle_params"]
C -->|"get_bminmax<br/>for all num_surf"| D["bmin/bmax arrays"]
D --> E["write_output"]
E -->|"generates"| F["bminmax.dat<br/>101 rows"]
G["class_parts.dat<br/>6 columns"] --> H["plot_classification.py"]
F --> H
H -->|"bin in Python"| I["class_jpar.pdf<br/>class_ideal.pdf"]
File Changes1. examples/classification/postprocess_class.f90
|
Code Review by Qodo
1. Plot binning normalization bug
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
Comment on lines
+10
to
+22
| hs = 1.0 / ns | ||
| pmin = 0.0 | ||
| pmax = np.max(perp_inv) | ||
| hp = (pmax - pmin) / nperp | ||
|
|
||
| prompt = np.zeros((nperp, ns)) | ||
| regular = np.zeros((nperp, ns)) | ||
| stochastic = np.zeros((nperp, ns)) | ||
|
|
||
| for ipart in range(len(s)): | ||
| i = min(ns, max(1, int(np.ceil(s[ipart] / hs)))) - 1 | ||
| k = min(nperp, max(1, int(np.ceil(perp_inv[ipart] / hp)))) - 1 | ||
|
|
There was a problem hiding this comment.
1. Plot binning normalization bug 🐞 Bug ✓ Correctness
bin_classification() normalizes perp_inv by its sample max, but doplot_inner() overlays bmin/bmax boundary curves in a different normalization, so the heatmap bins and boundary curves are on inconsistent y-scales and the plot can be physically misleading. If all perp_inv values are identical (e.g., all 0), hp becomes 0 and the script divides by zero when computing bin indices.
Agent Prompt
### Issue description
`examples/classification/plot_classification.py` bins `perp_inv` using `pmax = np.max(perp_inv)`, but overlays curves computed from `bminmax.dat` that assume normalization by the global minimum B (`bmin_global`). This makes the plotted heatmap and boundary curves inconsistent, and can also crash when `pmax == pmin` (hp=0).
### Issue Context
The boundary curves use `bmin_global / bminmax[:, 1|2]`, which is consistent with plotting `perp_inv * bmin_global` (dimensionless, ideally in [0,1]). The binning should use the same quantity, and must guard against `hp == 0`.
### Fix Focus Areas
- examples/classification/plot_classification.py[5-31]
- examples/classification/plot_classification.py[33-56]
### Suggested approach
- Load `bminmax` first in `main()` and compute `bmin_global`.
- Define `perp_norm = perp_inv * bmin_global` and clip to `[0, 1]` (or set explicit bounds).
- Update `bin_classification()` to bin `perp_norm` in a fixed range `[0,1]` (so `hp = 1.0/nperp`), avoiding `pmax=np.max(...)`.
- Add a defensive guard so `hp` cannot be 0 (even if you keep data-driven pmax).
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
Comment on lines
771
to
775
| !$omp critical | ||
| bmod = compute_bmod(z(1:3)) | ||
| if (num_surf > 1) then | ||
| if (num_surf /= 1) then | ||
| call get_bminmax(z(1), bmin, bmax) | ||
| end if |
There was a problem hiding this comment.
2. Classifier skips bminmax for num_surf=0 🐞 Bug ✓ Correctness
compute_pitch_angle_params() now calls get_bminmax when num_surf/=1, but trace_orbit_with_classifiers() still only calls get_bminmax when num_surf>1. If num_surf remains 0 (e.g., startmode=2 reading a start.dat with varying radii), classification can use stale bmin/bmax and compute incorrect passing/trap parameters.
Agent Prompt
### Issue description
`src/simple_main.f90` updated the bmin/bmax lookup guard to `num_surf /= 1`, but `src/classification.f90` still uses `num_surf > 1` before calling `get_bminmax`. This can skip per-radius bmin/bmax updates when `num_surf==0`.
### Issue Context
Classifier computations (`passing`, `trap_par`, `perp_inv`) depend on accurate `bmin/bmax`. For `num_surf=0` (volume / mixed-radius starts), bmin/bmax should be obtained via `get_bminmax`.
### Fix Focus Areas
- src/simple_main.f90[762-777]
- src/classification.f90[148-155]
### Suggested approach
- Change the guard in `trace_orbit_with_classifiers` from `if(num_surf > 1)` to `if (num_surf /= 1)` (matching simple_main), or explicitly handle `num_surf==0` as a case that must call `get_bminmax`.
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
krystophny
changed the title
krystophny
force-pushed
the
feature/venv-setup
branch
from
31cfad8 to
50a0fce
Compare
The error stop treating fast_class and tcut as mutually exclusive was wrong: fast_class enables J_parallel + topological classifiers while tcut enables the fractal dimension classifier. They are complementary. Also fix the classification example referencing QA instead of QH VMEC file (mismatched with the Makefile download target), and add a combined classifier golden record test exercising both features together.
The samplers-via-startmodes refactor (08cf86f, March 2025) broke classification diagnostics by removing volume sampling and bminmax output. This made it impossible to reproduce Figure 12(a-b) from the 2020 JPP paper, as reported by Rohan Ramasamy. Fixes: - startmode=1 with num_surf=0 now triggers volume sampling (uniform random s in [0,1]) instead of single-surface field line sampling - compute_pitch_angle_params calls get_bminmax for num_surf=0 (not just num_surf>1), giving per-particle bmin/bmax at varying radii - write_output generates bminmax.dat (101 rows) when classification is active, needed by the plot script - plot_classification.py reads class_parts.dat directly (all 6 columns) and bins particles in Python, replacing the standalone Fortran postprocessor that read only 5 of 6 columns - Fixed hardcoded bminmax[250,1] index that assumed 1001 rows - Deleted postprocess_class.f90, simplified Makefile
The binned array has shape (nperp, ns) which imshow displays correctly as y=J_perp, x=s. The erroneous .T transposed this to (ns, nperp), swapping the axes. Verified pixel-identical output against the original Fortran postprocessor.
Add single-surface (s=0.3, s=0.6) and volume input files for the QH configuration, a run script, and a plotting script that generates Figure 6 style loss-time vs trapping-parameter plots, Figure 8 style orbit classification in (theta, v_par/v) space, volume classification in (s, J_perp) space, and Table 1 style regular-fraction statistics. Also fix plot_orbit_chartmap_comparison.py to handle missing netCDF4 gracefully instead of crashing the test.
…example The fast_class feature incorrectly propagated the J_parallel classifier completion status (ierr_cot) to the integration error flag (ierr), causing ALL classified orbits to be marked as lost. Regular orbits classified by J_parallel should stop tracing early and be counted as confined; only chaotic orbits should continue tracing to the end. This produced wildly wrong confined fractions (e.g., fc=0.57 instead of fc=0.94 for QI s=0.3). The fix sets regular=.True. for J_parallel-classified regular orbits (when class_plot is off), which correctly: (a) stops integration early, (b) counts the orbit in confpart at all remaining time steps, and (c) sets times_lost to trace_time (confined). When class_plot is on, the old behavior is preserved so classification output is written at ntcut. The classification example is extended to reproduce all results from the JPP 2020 paper (Albert, Kasilov, Kernbichler) using the actual paper equilibria: QI (Drevlak 2014), QH (Drevlak 2018), QA (Henneberg 2019). Input configs for all cases (Figures 5-8, Table 1, volume classification) are provided under configs/. The plot script generates publication-quality figures matching the paper. Confined fractions validated against the paper at 1000 particles: QI s=0.3 fc=0.938, QH s=0.3 fc=0.993, QA s=0.6 fc=0.697.
Include the three VMEC equilibrium files used in the JPP 2020 paper: - QI (Drevlak 2014): wout_23_1900_fix_bdry.nc - QH (Drevlak 2018): wout_qh_8_7.nc - QA (Henneberg 2019): wout_henneberg_qa.nc Include reference output plots (100k particles, trace_time=1s): - fig5_losses_qi.pdf, fig6_losses_qh.pdf, fig7_losses_qa.pdf - fig8_classification.pdf, volume_classification.pdf Update run_all.sh to use repo-local data/ directory instead of hardcoded paths, so Rohan (or anyone) can reproduce by running: make && ./examples/classification/run_all.sh
- Remove wout .nc files and output PDFs from SIMPLE repo (belong in the private proxima-simple-classification repo) - run_all.sh now auto-downloads wout files from the proxima repo - Add qi_volume.in and qa_volume.in for volume classification of all 3 configs - Fix Fig 8 plot to use dense colored grid (imshow) matching the paper's "brazilian flag" style instead of sparse scatter markers - Update plot_paper_results.py to generate volume plots for all configs
Replace the sparse 100x100 dominant-category grid with a 50x50 Gaussian-smoothed fractional composition approach that fills the entire (theta/pi, v_par/v) plane with continuous colors -- the "brazilian flag" look matching the JPP 2020 paper Figure 8. Reduce overlay marker sizes for early-loss and chaotic particles.
Change num_surf > 1 to num_surf /= 1 in trace_orbit_with_classifiers, matching the fix already applied to compute_pitch_angle_params. For volume sampling (num_surf=0), each particle now gets bmin/bmax from its actual flux surface rather than using fixed values from sbeg.
Use 4 categories (passing/regular-trapped/chaotic/prompt-loss) instead of 3 so the trapped-passing boundary (dashed line) is visible as the gray-to-colored transition. Previously passing and regular-trapped were both yellow, hiding the boundary inside a uniform region.
Extract init_bminmax_arrays from get_bminmax lazy initialization into an explicit subroutine. Call it during startup while magfie is still in VMEC mode, ensuring find_bminmax scans geometric (theta, phi) angles rather than canonical coordinates.
Skip passing particles in the binning so the passing region is white. The trapped-passing boundary (dashed line) now visibly traces the edge of the colored data region.
krystophny
force-pushed
the
fix/classification-regression
branch
from
a3210dd to
cec6270
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Fixes multiple classification regressions and extends the classification example to reproduce all results from the JPP 2020 paper (Albert, Kasilov, Kernbichler).
Bug fix: fast_class marking regular orbits as lost (579d8ef)
The
fast_classfeature propagated the J_parallel classifier completion status (ierr_cot) directly to the integration error flag (ierr), causing ALL classified orbits -- including regular ones -- to be marked as lost. This produced wildly wrong confined fractions (e.g., fc=0.57 instead of fc=0.94 for QI s=0.3 with 100k particles).The fix sets
regular=.True.for J_parallel-classified regular orbits, which correctly: (a) stops integration early, (b) counts the orbit in confpart at all remaining time steps, and (c) recordstimes_lost = trace_time(confined). Whenclass_plotis on, the old termination-at-ntcut behavior is preserved so classification output is written.Earlier fixes in this branch:
Extended classification example reproduces Figures 5-8 and Table 1 using the actual paper equilibria:
Verification
Test passes after fix
Confined fractions match paper (1000 particles, same as paper)
Before fix (fast_class bug, 100k particles)
All trapped particles were marked as lost, giving fc ~ 0.57 (only passing particles counted as confined).
Generated plots
Figure 5 (QI): fig5_losses_qi.pdf
Figure 6 (QH): fig6_losses_qh.pdf
Figure 7 (QA): fig7_losses_qa.pdf
Figure 8 (classification): fig8_classification.pdf
Volume classification: volume_classification.pdf
Test plan
make test-fastpasses 100% (46/46)